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ABSTRACT 

As part of an ongoing research program at the Naval 
Postgraduate School concerned with aircraft vulnerability 
and, in particular, the "Hydraulic Ram" phenomenon, two 
Fortran IV computer codes were written to analyze in two 
dimensions problems concerning structure-fluid interaction. 
This thesis presents the method selected for solution, 
details the codes written and presents test cases and 
example problems. User's instructions and program listings 
for both codes are also included. 
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I. INTRODUCTION 



Many important structural problems in aerospace engi- 
neering Involve components that are in contact with a fluid. 
One example is the fuel tank located in a wing or fuselage 
of an aircraft. Another example is the fuel tank in a 
liquid fueled rocket. An example in the field of civil 
engineering is the dam and reservoir system. 

An ongoing research program at the Naval Postgraduate 
School is concerned with aircraft vulnerability in the fuel 
tank area. This program is currently devoted to a study of 
the Hydraulic Ram Problem. "Hydraulic Ram" is the dynamic 
loading of fuel tanks when impacted by bullets or warhead 
fragments. This loading gives the appearance of a simple 
pressure pulse with a resulting deformation or rupture of 
the tank. It has been found, however, that the loading and 
deformation are a combination of a number of events called 
"hydraulic ram components." [Ref. 1] 

As part of the theoretical analysis of the fuel tank 
response to a penetrating projectile two computer codes were 
written to solve for the response of the fluid and the walls 
of a idealized two-dimensional fuel tank. The codes will 
compute the pressure throughout the fluid and will calculate 
the wall response to this pressure. The user can prescribe 
the problem to be analyzed, including ) but not limited to: 
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Fluid surface pressure disturbances, initial wall displace- 
ments, initial wall velocities, pressure distributions through 
the fluid, etc. 

This thesis presents the method selected for the solution 
of problems involving two-dimensional fluid-structure inter- 
action. The types of problems which can be solved using the 
codes are described. Section II contains the theoretical 
description of the problem, and the following two sections 
outline two finite difference formulations that were used 
to obtain solutions. Section V contains a description of 
each code and describes the modifications necessary to 
change the boundary and initial conditions. Section VI 
presents the solutions obtained with the codes, including 
the check cases and several example programs, and discusses 
the comparison of results obtained from the two code formula- 
tions. Users’ instructions and program listings for both 
formulations are included in Appendices A and B. 
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II. DERIVATION OF GOVERNING EQUATIONS 



A. THE WAVE EQUATION FOR AN INVISCID FLUID 
1 . Equation of Motion 

Euler’s equation of motion can be given in the form 
P H + grad p = 0 (1) 



where P is the fluid density, p is the fluid hydrostatic 
pressure (positive in compression), q is the velocity vector 
and t is time. The total derivative — is given by 



D _ 9 , 9 t 9 , 9 

Dt ~ 3t U 3x V 3y W 3t 



( 2 ) 



where u, v, w are the components of the velocity in the 
x, y, z directions respectively. When the fluid velocity is 
sufficiently small, 



D_ 9_ 
Dt ~ 3 1 



( 3 ) 



Superscript dots will be used to denote partial derivatives 
with respect to time. 

2 . Constitutive Equation 

The constitutive equation of the inviscid fluid is 

given by 

p = - K div q (*0 

where K is the bulk modulus . 
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3 . Equation of Continuity 



p 5t + div q = 0 

Applying Eq . (3) we obtain 



( 5 ) 



— p + div q 
P p 



0 



( 6 ) 



Eliminating div q 



between Eqs . 



(4) and (6) gives 




(7) 



4 . Potential Equation 
Define a potential 

q - grad <j> (8) 

Thus Eq . (1) becomes 



p grad $ + grad p = 0 (9) 

when Eq . (3) is applied. Assuming p to be essentially 
constant throughout the fluid, Eq . (9) becomes 

grad (p4> + p) = 0 (10) 



or 



P4> + P 



constant 



( 11 ) 



The constant term is set to zero since it implies a base 
pressure which is constant throughout the domain. 

Differentiating Eq. (1) with respect to time gives 
• • 

pq + grad p = 0 (12) 

Applying Eqs. (*1) and (8) 

p grad <j> + grad (~K div grad <f>) = 0 (13) 

Since p is assumed to be essentially constant, Eq . (13) 
becomes 



grad(p<}> - K div grad <f>) = 0 



or 

.. p 

grad (pi}) — KV cf> ) ■ 0 



( 14 ) 



(15) 



Then 



or 



.. p 

p<j> - KV <f> = constant 



p4> = KV 2 <p 



( 16 ) 



(17) 



Recalling that the speed of sound in a fluid can be deter- 
mined from the relationship 



c 



2 



K 

P 



(18) 
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we obtain the wave equation 



c 2 V 2 (J> = <f> (19) 

B. THE GOVERNING EQUATION FOR THE STRUCTURE 
1 . Two-Dimensional Idealization 

Consider a long open box composed of five uniform 
plates. It can be idealized in two dimensions as a structure 
composed of three unit width beams attached at right angles 
as shown in Fig. 1. In the analysis each beam is considered 
separately, with its own set of independent boundary conditions. 
In all of the examples considered, each beam was assumed to 
be clamped at both ends . ' 

The equation for the' lateral displacement of a uniform 
beam under dynamic load is [Ref. 2] 

a 4 

El - — jj- w(x,t)+mw(x,t)=-p (20) 

9x 

where E is Young’s modulus, I is the section moment of 
inertia, m is the mass per unit length, w(x,t) is the lateral 
displacement of the beam, and p is the pressure applied to 
the beam. A positive pressure is in the opposite direction 
as w(x,t) as shown in Fig. 1 for the tank. 

C. INTERFACE CONDITIONS 

From Eq . (11) the pressure on the beam is the fluid 
pressure given by 
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P = - p4> 



( 21 ) 



where again a positive pressure is in the opposite direction 
as w(x,t) . 

Another interface condition is the requirement that the 
normal velocity of the fluid potential at the interface is 
equal to the velocity of the wall, i.e. 

w ( x , t ) = K? <J> (22) 

3n 

where — indicates a partial derivative normal to the wall. 
3n 

Where the fluid has a free surface and the pressure is 
zero, the boundary condition is given by 

<f> = 0 (23) 

When an excitation pressure is applied to the free surface 

* - - ( 2*0 



D. SUMMARY 

In summary, the equations used in the analysis are 



c 2 V 2 <{> = <J> 



(19) 



El 



P = 



w(x 



3 

— jr w(x,t) + mw(x,t) = -p (20) 

3x 




( 21 ) 




3n 



( 22 ) 
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III. NUMERICAL METHOD OP ANALYSIS - ACCELERATION FORMULATION 



A. FINITE DIFFERENCING SCHEME 
1. The Fluid 

The wave equation is hyperbolic in nature and for 
this reason an explicit formulation for the timewise response 
is chosen. The conventional central finite difference 
equations are used in both space and time, and when applied 
to Eq. (19) give [Ref. 3] 



(Ax)' 



{♦i.j- 



1 + *1,0+1 4 *i,j + *i+l , j + * 



i-1 3 0 } 



( 5F )2 K-l - 2 *k + *k + i} i,J 



(25) 



where Ax = Ay in a square mesh, i and j are mesh point 
reference indices, and k is the time step. <)k ^ ^ + - L is the 
only unknown. All other terms are known or can be evaluated 
prior to the calculation of <j>^ j . 

2 . The Beam Equation 

Conventional central finite difference schemes are 
also used for the walls in both space and time. The resulting 
equation is given by [Ref. 3] 
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( 26 ) 



El 

(Ax) 1 * 




^ w i_l + ^ w i 





+ 



m 

(At) 2 




2w, + w. 



k+1 



-P 



i,k 



where 

Pi.k - - 2St {-Vi + w}, . (27) 

and j = 1 or j = J max . See Fig. 1. 

A similar pair of equations can be written for the wall 

along the boundary where i = 1. In Eq. (26) w. , , , is the 

1 y Kt J. 

only unknown in the explicit scheme. Note that the pressure 
given by Eq. (27) is a function of the fluid potential at 
the wall . 

3 • The Calculation Scheme 

When the two systems are combined, the other linking 
mechanism is the interface condition on the fluid and wall 
velocities. Using central finite differences, Eq . (22) is 
given by [Ref. 3] 

2At { -w i , k-1 + w i,k+l} = 2A7 {“^i,.]-! ,k + *i,j+l,k} (28) 
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In the calculation scheme the wall deflection at each new 



step k+1 is computed first, and then the velocity potential 
is computed throughout the fluid. In order to 'compute the 
wall displacement at k+1 the pressure must be known at the 
preceding time step k, which requires knowledge of the 
velocity potential at the wall at the new time step, k+1, 
according to Eq . (27). The solution to this problem is to 
solve for the wall displacement by implicitly solving for 
the new velocity potential at the wall. Thus, using 
Eq. (25) to eliminate <f>. . , ,, from Eq . (27) > and using 
Eq. (28) to eliminate <b. . , . (when j = l ) from the result 

1 % J — X 3 K 



A similar equation can be obtained for the wall where 

j = j . Substituting p given by Eq . (29) into Eq. (26), 
max 

combining terms, and solving for w^ gives 
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w 



i,k+l 



25 ht K.j.k - 2<, l ; j,k-l + 



/ cAt v 2 
^ Ax ; 



[<f, i-l,j ,k + *i+l,l,k 



4< ^i,l 3 k + 24> i,2,k 



+ 



Ax 

At 



w i,k 



4 



+ 




w. 






El 

m 



(- A — 9 ) 2 (w. „ . - 4w. n . + 6w. - w. ,, 

(Ax) 2 l 1-2 > k i-l, k i,k l+l, 

'i+2,k} 



k 



+ w. 



/ t i 1 At 2^ 

/ (1 + — 7 C ) 

m Ax 



(30) 



Equation (30) is written for the wall where j = 1 , 0 <_ x < aH. 
For x > aH, p = 0. Once the wall deflection is known at k+1, 
the interface conditions on the fluid are known. The calcula- 
tion of the fluid velocity potential is now quite straight- 
forward using Eq . (25) and Eq . (28) at an interface. 

On the free surface of the fluid the velocity 
potential may be excited by an applied pressure. In this 
case, applying central finite differences to Eq . (2*0 gives 



w ■ Vi - 2 f <p b 5 k (31) 

where (P^)^ is the pressure exerted on the boundary at time 
step k. 
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B. INITIAL CONDITIONS 



1. The Fluid 

Since the governing partial differential equation 
(Eq. (19)) is second order with respect to time, two initial 
conditions are required; 1) the velocity potential at time 
zero', and 2) the first time derivative of the velocity 
potential, which is proportional to any initial pressure in 
the fluid. 

2. The Wall 

The beam equation (Eq. (20)) is also second order 
with respect to time, and two initial conditions are required; 
1) the initial displacement of the wall, and 2) the initial 
velocity of the wall. 

C. THE TIME STEP FOR NUMERICAL STABILITY 
1. The Fluid 

For hyperbolic systems solved using Eq . (25), the 
non-dimensional time step ratio 



r = 



At 1 

Ax' 



(32) 



where primes denote non-dimensional quantities, must be less 
than 1 for numerical stability for a one-dimensional system 
[Ref. 3 ] and must be less than \T 2^ /2 for numerical 
stability for a two-dimensional system [Ref. 4]. 



20 



Non-dimensional time t* is given by 



t» 



tc 

L 



( 33 ) 



where L is a, reference length. The non-dimensional length 
x* is given by 




(3*0 



Substituting Eqs . (33) and (3**) into Eq . (32) gives 



t < 



Ax 



JT 



(35) 



for numerical stability of the two-dimensional system. 

2. The Wall 

Non-dimensionalizing the time and length coordinates 
in the beam equation gives [Ref. 3] 



t ' » 4 




(36) 



x' 



x 

L 



(37) 



The beam equation is parabolic in nature [Ref. 3] and 
requires a non-dimensional time step ratio 



At * 1 

(Ax') 2 " 2 



(38) 
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for numerical stability. Substituting Eqs. ( 36 ) and (37) 
into Eq. ( 38 ) gives 




(39) 



Therefore 



At 




(40) 



3 . The Combined System 

The assumption has been made that the fluid will 
require the smallest time step for stability. Thus, there 
is a restriction on the value of the parameters used in the 
combined analysis, in particular the physical parameters of 
the beam. This restriction is given by 




(41) 



according to Eqs. (35) and (40). 
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IV. NUMERICAL METHOD OF ANALYSIS - VELOCITY FORMULATION 



A. 'VELOCITY FORMULATION APPROACH 

Kreig and Monteith [Ref. 5] have suggested a modification 
to the timewise finite differencing scheme outlined in 
Section III. This modification is described as a "velocity" 
formulation . 

Consider their example of the one-dimensional spring mass 
system. 



Mw + Kw 



0 



(42) 



Using central finite differences 3 Eq . (42) becomes 



(At) 2 



M 



{ 




+ Kw. 



k 



0 



(43) 



or 




(44) 



Now consider the velocity formulation 



Mv + Kw = 0 



w 



V 



(45) 
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Defining the velocity at half time steps and the displacement 
at whole time steps leads to 



At 



f 



k-h 



+ v k + %} 



KW k = 0 



(H6) 



or 



v k +% = v k- !s - 1 At w k 



(*»7) 



and 



w 



k+1 



= w, + At v. 



k+% 



(48) 



With either formulation the time step is the same, and in 
the real number system without round-off error the results 
would be identical. However, computers do not compute in 
the real number system. Round-off errors are always present. 
Consider the term in brackets in Eq . (44) 



[2 - | (At) 2 ] 



Depending on the number of digits retained, the significance 
K p 

of jjj (At) could be completely lost in round off. 



B. FINITE DIFFERENCING SCHEME 
1 . The Fluid 

The fluid equation derived in Section III can be 
written as 



2 2 • 
c v <J> » ip 



(*»9) 



<f> = 'I' 



(50) 



A central finite difference scheme is proposed centering 
the "velocity" at whole time steps and the "displacement" 
at half time steps. Note, the pressure is directly computed 
at each whole time step since 



" P 



(51) 



The finite difference equation for the fluid is given by 



(A 



x ) 2 + j + *i+l,J + ‘t'i-l. j} k _ Ss 



■& {-*k-l + *k} . 

1 i « 



(52) 



where \p. . . is the only unknown. The potential is obtained 

1 , J j K 

from the finite difference expression for Eq. (50) 



*l,J,k+is ~^i j J >k-k 



At *i,J,k 



(53) 
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2. The Wall 



Similarly, the beam equation can be written 

El — Jr w(x,t)+mv = -p (54) 

3x 

and 

w(x,t) = v (55) 



For these equations a central finite difference scheme is 
proposed centering the velocity v at half time steps and 
the displacement w at whole time steps. Note that the 
velocity needed for the interface condition is directly 
computed. 



i,k-*s Ax 



{'♦i.J-i + hj+iL, 



(56) 



The finite difference equation for the beam is 



— ^ { Wl _ 2 - Hw i _ 1 + 6w ± - 4w 1+1 + w 1+2 }^ + 

SF + V i,k + J = 

The unknown is v. , ,, . From Eq . (55) 

1 y rC > ^ 

w i,k+l = w i,k + At v i,k+*s 



(57) 



(58) 
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3 . The Combined System 



With the velocity function there is no need to solve 
for the beam displacement in terms of an implicit solution 
of the potential at the wall. In the velocity formulation 
the pressure is computed directly; no differencing is 
required. Figure 2 graphically illustrates the time step 
centering of both the acceleration formulation and the 
velocity formulation. 

The solution scheme for the velocity formulation is: 

1. Solve for the wall velocity, Eq . (57)> based on the 
known pressure at the wall and the previous displace- 
ment and velocity. 

2. Calculate the new wall displacement, Eq . (58), based 

on the new wall velocity and the previous displacement. 

3. Solve for \p, Eq. (52), based on the previous the 
new wall velocity (for the interface conditions) and 
the previous ([>. 

4. Calculate the new 4>, Eq . (53) , based on the previous 
(f> and the newly calculated \p . 

5. Return to step 1 and repeat. 

As in the example given by Kreig and Monteith this formula- 
tion is exactly equal to the acceleration formulation if the 
beam and fluid equations are considered separately. However, 
in the combined system the fluid equation is displaced 
exactly one-half time step from the fluid equation used in 
the acceleration formulation. In the acceleration formulation 
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ACCELERATION FORMULATION 



<f> <f> 

W W 

> TIME 

-1 0 +1 +2 + 3 . 

I T I 

£ 

W 

VELOCITY FORMULATION 




FIG. 2 



28 



the pressure depends on (j^+i an ^ 4 ) i c _x whereas in the velocity 
formulation, it depends on and 4> lc+ ^ • Thus, even with 

no round-off error there will be a difference between the two 
solutions. Refer to Part D for a discussion of the signif- 
icance of this difference. 



C. INITIAL CONDITIONS 

As with the acceleration formulation a total of four 

initial conditions are required. However, recall that two 

of the four variables are defined one-half time step on 

either side of t = 0. They are the wall velocity and the 

velocity potential of <}> in the fluid. 

The actual initial conditions can be denoted as w and 

o 

<j> 0 . The computer code assumes that 



v 




w,. 



( 59 ) 



and 



<f> 



*o 



( 60 ) 



D. THE TIME STEP FOR NUMERICAL STABILITY 

Theoretically, since the same parameters are used in 
both formulations, the time step for numerical stability 
should be exactly the same for both formulations. In fact, 
this was found to be true for the beam equation and for the 
fluid equation when each was tested separately. However, 
when the systems were coupled, and. the maximum theoretical 



29 



time step was used, substantial instability exhibited 
itself immediately in the computer results. The instability 
occurred at the interface and it occurred within two time 
steps of a pressure pulse striking a wall. A subsequent 
velocity was imparted to the wall which was always large 
enough to cause a reversal in the fluid pressure at the wall 
at the next time step. 

A primitive but effective solution to the instability 
problem was obtained by reducing the time step arbitrarily 
to 



At 



Ax 

c ry 



(6l) 
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V. DESCRIPTION OF THE COMPUTER PROGRAMS 



Two Fortran IV programs were written and tested. The 
theory behind each is outlined in Sections III and IV above. 

The two programs are quite similar in construction and 
require similar amounts of core storage. C.P.U. time 
required is also similar. Only one program will be described 
in detail, but differences between the two will be pointed 
out where they occur. 

A. GENERAL FEATURES 

The programs allow for complete flexibility in any of 
the possible initial conditions, both on the walls and on 
the fluid. At present the program allows for three walls 
and a variable fluid depth. The ends of each wall are 
considered to be clamped. The wall boundary conditions can 
be changed by any future user by changing six equations in 
the velocity formulation and twelve in the acceleration 
formulation. The location of these equations will be pointed 
out below. The program provides for a non-dimensional thirty 
by thirty mesh point domain, with the reference length being 
the length of the horizontal wall (See Fig. 1). All variables, 
vectors and arrays are placed in blank Common, to allow for 
complete flexibility in communication between subroutines. 
Further all vectors, and arrays are initialized to zero value 
before analysis begins. 
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Subroutines are provided to: 

1. Compute constants derived from input parameters, 
normalize the dimensions of the problem, and compute 
or set initial conditions. 

2. Output data obtained for each mesh point in the domain. 
The format is tabular. 

3. Compute the beam equations. 

4. Compute the fluid equations. 

5. Provide for a special graphical output. 

Each subroutine will be described in detail. 



B. 



MAIN PROGRAM 
1 . Variable 
WA (I ) 

WB(I) 

WC(I) 

VA(I) 

VB(I) 

VC (I ) 

UA (I ) 

UB( I) 

UC (I ) 

XA(I) 

XB(I) 

XC(I) 

WADOT(I) 

WBDOT ( I ) 



Description — Velocity Formulation 
left vertical wall displacement at k 
horizontal wall displacement at k 
right vertical wall displacement at k 
left vertical wall displacement at k+1 
horizontal wall displacement at k+1 
right vertical wall displacement at k+1 
left vertical wall velocity at k-h 
horizontal wall velocity at k-h 
right vertical wall velocity at k-h 
left vertical wall velocity at k+h 
horizontal wall velocity at k+h 
right vertical wall velocity at k+h 
initial left vertical wall velocity 
initial horizontal wall velocity 
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WCDOT(I) 


- initial right vertical wall velocity 


PRA ( I ) 


- initial left vertical wall pressure 


PRB(I) 


- initial horizontal wall pressure 


PRC (I ) 


- initial right vertical wall pressure 


PK (I ) 


- fluid free surface applied pressure 


IWA(I) 


- integer value of left vertical wall 
displacement 


IWB(I) 


- integer value of horizontal wall 
displacement 


IWC(I) 


- integer value of right vertical wall 
displacement 


P(I,J) 


- value of 4> in the fluid at k-h 


PN(I,J) 


- value of (J) in the fluid at k+% 


PL(I,J) 


- value of \p in the fluid at k-1 


PRESS ( I, J) 


- value of pressure in fluid at k 


PHIDOT(I,J) 


- initial value of ip in the fluid 


IPRESS(I ,J) 


- integer value of pressure in the fluid. 



In the displacement formulation the following variables differ 



va,vb,vc 


are not used 


WA,WB,WC 


are values of wall displacement at k 


XA,XB,XC 


are values of wall displacement at k+1 


UA,UB,UC 


are values of wall displacement at k-1 


P(I,J) 


are values of $ at k 


PN(I,J) 


are values of $ at k+1 


PL(I, J) 


are values of <j> at k-1 
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2 . Main Program Functions 



In the velocity formulation the main program 
accomplishes the following functions: 

1. Controls calls to other subroutines. 

2. Computes first value of PN(I,J) from initial conditions. 

3. Saves the value of the displacement at a selected 
wall mesh point at selected time steps for later 
plotting as a time vs. displacement graph. 

*1. Turns off a free surface pressure disturbance at any 
preselected time step. 

In addition the "acceleration" formulation program 

5. Computes the k = 1 value of wall displacement from 
initial conditions. 

C. SUBROUTINE INP(IER) 

This subroutine 

1. Reads all input data. 

2. Computes constants. 

3. Assigns values to initial conditions other than zero 
value . 

1 . Variable Description 



RHOF 


- Density of the fluid 


E 


- Young's Modulus 


RL 


- Horizontal length 


RH 


- Vertical height where RH <_ RL 


AL 


- Fractional depth of fluid 0 < AL 5 1.0. 


C 


- Speed of sound in the fluid. 






TH 


- Thickness of beam 


RHO 


- Density of beam material in lbs 
per ind 


NH1 


- Number of vertical mesh points 


N1 


- Number of horizontal mesh points 
where NH1 £ N1 = 30 


DELX 


- Width of mesh 


DELT 


- Time step 


RIO 


- Computed value of I, moment of 
inertia 


NA 


- Number of vertical mesh points in 
the fluid . 



C1,C2,C3 - Constants common to both programs 

C4 , C5 j C 6 , C7 , C 8 , C9 j D - Constants used only in acceleration 

formulation 

D1,D2,D3 - Constants used In velocity formulation 



NH 


- Equals NHI - 1 


NHI 


- Equals NHI - 2 


N 


- Equals N 1 - 1 


NAI 


- Equals NA - 1 


2 . 


Input Cards 

Only two input cards are used to input RL, RH, AL, 



C, TH, RHO, RHOF, E. The format will be described in 
Appendix A. 

3 . Non-Zero Initial Conditions 

Any non-zero initial conditions are inserted by the 
user at the end of Subroutine INP . 
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D. SUBROUTINE OUTP(IER) 

Subroutine OUTP prints out in tabular format the values 
of each wall mesh point displacement and the pressure at 
each mesh point throughout the domain. 

E. SUBROUTINE WALL(IER) 

Subroutine WALL computes the values of the wall velocity 
and displacement in the velocity formulation and the wall 
displacement in the acceleration formulation. If a future 
user desired to change the boundary conditions of the beam 
equation, then any equation with a computed answer stored 
in a variable indexed with a 2, N, or NH would have to be 
changed. In the velocity formulation only the equations 
for VA ( 2 ) , VA(NH), VB(2), VB(N) , VC(2), VC(NH) need to 
be changed. In the acceleration formulation these equations 
are for XA(2), XA(NH), XB(2), XB(N), XC(2), XC(NH) and recall 
that these equations must also be changed in the main program. 

F. SUBROUTINE FLUID(IER) 

Subroutine FLUID computes the value of rp and <f> at k and 
k+h respectively in the velocity formulation and computes 
the values of 4> at k+1 in the acceleration formulation. 
Pressure at each mesh point is also computed in FLUID. 

G. SUBROUTINE GRAPHO (IER) 

Subroutine GRAPHO is an attempt to make some sense out 
of the possible 1900 numerical values calculated at each 
time step. The subroutine locates an absolute maximum value 
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for wall displacement and fluid pressure and from these 
values normalizes the entire domain to integer values from 
-99 to +100. This information is then presented in a format 
quite similar to Fig. (1) . An integer value is placed at 
each mesh point in the domain. 



37 



VI. SOLUTIONS 



A. CHECK CASES 

In order to determine if the programs correctly modeled 
the system described above, several test cases were run on 
each formulation. 

1. The Fluid 

By controlling the subroutines called by the main 
program, the user is able selectively to compute results 
for the fluid alone or the wall alone. One test case for 
the fluid alone was a step pressure disturbance applied 
uniformly over the entire upper surface for a finite time 
interval. The walls were treated as rigid. This test 
showed that the pulse propagation and reflection behaved 
as expected in both formulations. The uniform pulse traveled 
uniformly through the fluid. The value of pressure doubled 
as expected at the horizontal rigid wall and returned 
uniformly to the free surface. See Pigs. 3 and 4. 

Another test case consisting of a point pressure 
disturbance at the middle of the free surface showed that 
two-dimensional pulse propagation in the fluid was as 

i 

expected. The amplitude of the pressure decreased radially 
from the point source. 

In both test cases the theoretical time step for 
stability proved to be valid in both formulations. 



38 



-o- THEORETICAL VALUE 
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-O- ACCELERATION FORMULATION 



PRESSURE 
IN PSI 



5 

-4- 



10 

-+- 



15 



20 

-4- 



25 

I ► 



<i: 

cd 



a 



QJ 



- 0 < 1 - 



FLUID 

SURFACE 



1 



<t> 



d 



□ < 

<i □ 



HORIZONTAL 



o-<&- 



TIME= .2167*10-3 S ec 
FIG. 5 



WALL 



39 



O- THEORETICAL VALUE 
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2. The Wall 



A free vibration analysis was conducted using the 
beam equation to test the wall response alone. The wall 
was initially displaced to conform to the first mode of 
vibration. Rogers [Ref. 1] gives a value of 51.3 cycles per 
second for the value of the frequency of the first mode of 
vibration for the beam parameters used in the test. The 
results from both formulations returned a first mode of vibra- 
tion frequency of approximately 50 cycles per second. Again 
the theoretical time step for stability proved to be valid. 

B. EXAMPLE PROBLEMS 

The types of problems these two codes were developed to 
solve involve the combined system. The first example prob- 
lem consisted of a uniform pressure pulse at the free sur- 
face of a flexible tank for 2.26 milleseconds . When a solu- 
tion was first attempted using the velocity formulation, the 
previously mentioned instability was discovered when the 
theoretical time step for numerical stability of the fluid 
was used. Since no other solutions to this type of problem 
were available, first priority centered on determining the 
cause and cure of the instability. A reduction in the time 
step as discussed in Section IV solved the instability prob- 
lem. The velocity formulation was run with a reduced time 
step, and the acceleration formulation was run using the 
theoretical time step. The results for the displacement 
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of the center of the horizontal wall are shown in Pig. (5). 
The dominant frequency of vibration is approximately 12 
cycles per second for both formulations. There is a small 
amplitude discrepency of about 5 % between the two formula- 
tions. The resulting vibration appeared to involve all 
modes in the beam. 

The acceleration formulation was also used to solve the 

same problem with a uniform pulse held for .226 milleseconds 
t h 

or 1/lCr the previous example. Figure (5) again denotes 
a frequency of approximately 12 cycles per second, and the 
amplitude is proportional to the previous case . Maximum 
amplitude is about one-tenth that of the previous case. 

It is interesting to note that if the fluid is considered 
as incompressible, and its mass is added to the mass of the 
beam, the frequency of free vibration as given by Rogers 
[Ref. 2] is 11.1 cycles per second. 

Another example problem was run with a pressure exerted 
at the two center mesh points of the free surface for approx- 
imately 450 microseconds. In this case, the first mode of 
vibration predominated in the horizontal wall at a frequency 
of 13 cycles per second. Again only minimal differences 
in amplitude were observed between the two formulations. 



DISPLACEMENT OP CENTER OP 
HORIZONTAL WALL IN INCHES 
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VII. CONCLUSIONS 



The results obtained to date appear to support the 
following views: 

1. The amplitude of vibration in the walls are depen- 
dent upon the length of time the free surface excitation 
pressure is applied. This would be expected purely from 
energy concepts. 

2. Choice of formulation is up to the user. Each has 
its own virtues and its vices. The velocity formulation 
requires less core storage and less CPU time for the same 
number of time steps. (The acceleration formulation requires 
15$ more CPU time than the velocity formulation; however, 

the time step required is 19 percent smaller in the velocity 
formuation.) The velocity formulation would be easier to 
modify . 

3. No claims are made as to the efficiency of either 
program as presently written. In fact it is estimated that 
the program size could be reduced considerably. 

4 . More work needs to be done in checking out other 
simple cases to prove the validity of the codes, for exam- 
ple, pressure pulses passing through the fluid at angles 
from the horizontal. 

5. Either program will give a future user valuable 
insight about structural behavior in contact with a fluid. 
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USER'S INSTRUCTIONS 

Input of initial conditions has already been discussed 
in Section V above. Actual input data is contained in two 
cards with the following format : 

Card 1 



Column 


Variable 


Definition 


Format 


Card 1 : 


1-10 


RL 


Reference length 


F10 .4 


11-20 


RH 


Vertical height 


F10.4 


21-30 


AL 


Fractional depth 
of fluid 


F10.4 


31-40 


C 


Velocity of sound 
in the fluid 


F10.4 


41-50 


TH 


Thickness of wall 


F10 .4 


51-60 


RHO 


Density of wall _ 
material in #/in^ 


F10.4 


Card 2: 


1-12 


RHOF 


Density of fluid 
in #-sec^/in4 


E12.4' 


13-24 


E 


Young's Modulus of 
wall material 


E12.4 
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